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ABSTRACT 

Bright circumstellar nebulae around massive stars are potentially useful to derive 
time-dependent mass-loss rates and hence constrain the evolution of the central 
stars. A key case in this context is the relatively young ejection-type nebula Ml-67 
around the runaway Population I Wolf-Rayet star WR124 (= 209 BAG), which 
exhibits a WN8 spectrum. With HST-WFPC2 we have obtained a deep, Ha image 
of Ml-67. This image shows a wealth of complex detail which was briefly presented 
previously by Grosdidier et al. (1998). With the interferometer of the Universite 
Laval (Quebec, Ganada), we have obtained complementary Fabry-Perot Ha data using 
GFHT MOS/SIS. 

From these data Ml-67 appears more-or-less as a spherical (or elliptical, with the 
major axis along the line of sight), thick, shell seen almost exactly along its direction 
of rapid spatial motion away from the observer in the ISM. However, a simple thick 
shell by itself would not explain the observed multiple radial velocities along the line of 
sight. This velocity dispersion leads one to consider Ml-67 as a thick accelerating shell. 
Given the extreme perturbations of the velocity field in Ml-67, it is virtually impossible 
to measure any systematic impact of the present WR (or previous LBV) wind on the 
nebular structure. The irregular nature of the velocity field is likely due to either large 
variations in the density distribution of the ambient ISM, or large variations in the 
central star mass-loss history. In addition, either from the density field or the velocity 
field, we find no clear evidence for a bipolar outfiow, as was claimed in other studies. 
On the deep Ha image we have performed continuous wavelet transforms to isolate 
stochastic structures of different characteristic size and look for scaling laws. Small- 
scale wavelet coefficients show that the density field of Ml-67 is remarkably structured 
in chaotically (or possibly radially) oriented filaments everywhere in the nebula. We 
draw attention to a short, marginally inertial range at the smallest scales (6.7-15.0 
xlO~^ pc), which can be attributed to turbulence in the nebula, and a strong scale 



break at larger scales. Examination of the structure functions for different orders 
shows that the turbulent regime may be intermittent. 

Using our Fabry- Perot interferograms, we also present an investigation of the 
statistical properties of fluctuating gas motions using structure functions traced by 
Ha emission-line centroid velocities. We find that there is a clear correlation at scales 
0.02-0.22 pc between the mean quadratic differences of radial velocities and distance 
over the surface of the nebula. This implies that the velocity field shows an inertial 
range likely related to turbulence, though not coincident with the small inertial range 
detected from the density field. The first and second order moments of the velocity 
increments are found to scale as {\Av(r)\) ~ r^'^ and {\Av{r)\'^) ~ r^'^. The former 
scaling law strongly suggests that supersonic, compressible turbulence is at play in the 
nebula, on the other hand, the latter scaling law agrees very well with Larson-type 
laws for velocity turbulence. Examination of the structure functions for different orders 
shows that the turbulent regime is slightly intermittent and highly multifractal with 
universal multifractal indexes a ~ 1.90-1.92 and Ci ~ 0.04 it 0.01. 

Subject headings: turbulence — ISM: bubbles — ISM: individual (Ml-67) — ISM: 
kinematics and dynamics — stars: individual (WR124) — stars: Wolf-Rayet 



1. Introduction 

1.1. Nebulae Ejected by Hot Stars 

Nebulae around central stars of planetary nebulae (CSPN) or Population I Wolf-Rayet (WR) 
stars have proven to be a powerful diagnostic tool to understand how mass is lost in the post 
main-sequence phases of stellar evolution. CSPN and some (~ 1/3; Marston 1995, 1996, 1999) 
massive WR stars are surrounded by circumstellar nebulae (CN) made up of stellar material 
ejected during successive stellar wind phases. On the whole, these nebulae often appear to 
have similar physical properties, which point to a common formation mechanism (Chu 1993) 
independent of the strong differences between both types of hot stars. 

The dynamics of planetary nebulae and CN surrounding Population I WR stars are extremely 
sensitive to the stellar wind velocity and mass-loss history. In the context of the interacting stellar 
wind model, a spherically symmetric, fast, hot stellar wind {vexp ^ 10^ km s""*^) is catching up 
and colliding with a pre-existing slower {vexp '^ 10"*^^^ km s^^), denser wind (CSPN: Kwok et al. 
1978; Chu et al. 1991; Balick 1994; Pop. I WR stars: Marston 1995, 1996, 1999; Garcia-Segura k 
Mac Low 1995; Garcia-Segura, Mac Low & Langer 1996, and references therein). The non-circular 
shapes of many CN suggest that the fast wind is sweeping up a slow wind that is flattened by 
rotation or binary effects (Garcia-Segura &: Mac Low 1995; Mellema & Frank 1995). However, a 
flattened, fast wind blowing out a spherically symmetric slow wind can also lead to a non-circular 
CN (Frank et al. 1998). 

Rayleigh- Taylor instabilities are expected to occur as a consequence of such interaction, as well as 
other hydrodynamical instabilities which may lead to turbulence in the nebula. On a qualitative 
basis, the interacting stellar wind model agrees with the observations, but it is still difficult to 
determine the precise initial, overall geometry of the colliding winds. In particular, the question 
arises as to the impact of pre-collision, apparently universal, clumping or possible turbulence of the 
hot-stai winds (Robert 1992; Balick et al. 1996; Lepine et al. 1996; Acker et al. 1997; Grosdidier 
et al. 1997; Eversberg et al. 1998; Lepine & Moffat 1999; Grosdidier et al. 2000, 2001; and 
references therein) on the dynamics/morphology of CN. Moreover, despite preliminary work done 



by Stone, Xu &: Mundy (1995), the incidence of a highly variable (as well in the spatial as in the 
time domain) mechanical luminosity [oc M{r^ 0, (j), t) x w?;o«i(^' ^' 't'l 0] ^till has to be investigated in 
detail. In particular, the turbulence of the ejected nebulae has to be studied in detail in order to 
test whether it is related to pre-existing turbulence (or any other process of fragmentation, hence 
variability) that already arose in the winds themselves before they interacted. 
Recall that the spectral variability originating in massive WR stars and CSPN shows the WR 
phenomenon to be remarkably similar in both (Grosdidier et al. 2000, 2001). Despite the strong 
differences between both types of hot stars, this supports the understanding of the WR spectral 
phenomenon as being purely atmospheric. This, in line with the believed similarity of the processes 
that are at play in the formation of the related CN (Chu 1993), places one in a position to check 
also for universality of the shaping of CN surrounding any hot stars. 



1.2. Turbulence in Ejected Nebulae and HII Regions 

Like HII regions, gas motions in nebulae ejected by hot stars are characterized by very large 
Reynolds numbers. For example, the Reynolds numbers exhibited by typical HII regions may 
be as large as 10^ (Boily 1993), using 1 pc for the characteristic spatial scale, 10 km s~^ for 
the velocity scale and considering the normal kinematic viscosity for ionized hydrogen gas at 
electronic temperatures of about 10,000 K. Such values are well above the critical values (10^-10^) 
determining the transition to turbulence. Like HII regions, it is expected that the kinematics of 
nebulae ejected by hot stars would be best described in terms of turbulent flows. Turbulence is 
characterized by a high degree of non-linearity in the governing equations and complex coupling 
mechanisms, leading to a variety of physical/transport processes occuring over a great variety of 
scale lengths. 

Turbulence is inherently apparently stochastic, since the flow variables {e.g. velocity, pressure, 
density) fluctuate in an unpredictable manner about their mean values. However, turbulence 
also exhibits specific length scales and scaling laws, with rapidly increasing numbers of smaller 
features as a result of energy cascading down a ladder leading ultimately to dissipation (Fleck 
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1996). At small scales, turbulent energy returns to the ambient medium as thermal energy. The 
nearly dissipationless energy cascade makes the flow not strictly unpredictable if one attempts to 
study quantitatively its structure via a statistical approach. Indeed, the signature of turbulence is 
confirmed by the presence of correlations in the flow variables between neighboring points that can 
be detected by means of statistical functions (Scalo 1984). The non-linearity of the equations of 
hydrodynamics causes high amplitude structures to change their shape in a way that continuously 
broadens the wave spectrum to include smaller and smaller scales. A simple, unique wave train is 
expected to steepen into a shock front damping at the leading edge. However, an intricate pattern 
of motions likely drives a cascade of energy dissipation taking on a fractal structure (Henriksen 
1994, and references therein). 

A statistical search for any degree of correlation between neighboring points could be carried out 
via two main approaches: autocorrelation functions or structure functions (these two statistical 
tools are more reliable than the dispersion-scale relation; Scalo 1984, Miville-Deschenes & Joncas 
1995). The so-called 'structure function' is defined as the average of the quadratic differences of 
velocities as a function of the spatial separation between the points where the velocity is measured 
(Miesch &: Bally 1994). In this study we favor the structure function calculation because it is 
more robust than the autocorrelation method. The robustness of the structure function technique 
originates in one important fact: its application requires the use of velocity differences rather 
than velocity 'coherence' (in the case of the autocorrelation method) between points separated 
by a given distance. Indeed, the use of the autocorrelation function assumes that the statistical 
properties of the fiows are not a function of position, the velocity and density fields being 
considered stationary. When using structure functions, we assume a less stringent hypothesis: 
local stationarity. However, note that at large scales the fiows always lose their stationarity (Scalo 
1984). The aim of the present paper is to describe the turbulent behavior of nebulae ejected by hot 
stars, both from the density structure (via high-resolution imagery) and the kinematical structure 
(via high spatial resolution Fabry- Perot interferometry) . The search for turbulence will be carried 
out through the calculation of the moments of the density /velocity increments for different orders 
(i.e. we will not limit our investigation to the second order structure function, as is generally 



the case in other studies) . Examination of the structure function for different orders additionally 
provides a way to estimate the level of intermittency in the flows. 

Ml-67, surrounding the WR star WR124 (WN8), was chosen for its relatively large apparent 
angular diameter (which allows one to study any possible energy cascade within a large span 
of scale lengths) and its suspected weak interaction with the ambient interstellar medium. The 
latter condition ensures that 1) the ejecta would be less affected by collision and mixing with the 
pre-existing ambient gas, and therefore 2) the dynamics and structure of the ejected nebula would 
best reflect the initial conditions of the ejected winds. This, in line with the apparently universal, 
ubiquitous clumping of hot stellar winds, places one in a position to test the possible impact of 
the wind fragmentation on the nebular fragmentation. Note that high resolution optical spectra 
of WR124 recently obtained by Acker (priv. comm.) at the Observatoire de Haute Provence 
(Prance) and Marchenko et al. (1998) results indeed show the central star of Ml-67 to be variable. 

2. The Nebula Ml-67 and its Central Star WR124 

Ground-based, narrowband images and discussion of the ejected nebulosity Ml-67 surrounding 
the cool nitrogen WN8 WR star WR124 (= 209 BAG) are given in Esteban et al. (1991), Esteban 
et al. (1993), Sirianni et al. (1998), Growther et al. (1999), and references therein. Until 1991 the 
central star of Ml-67 was often mistaken for a GSPN. From a study of the Na I D2 interstellar 
absorption line in the context of Galactic rotation, Grawford & Barlow (1991) showed that the 
distance to WR124 is 4-5 kpc. For such a large distance, the central star must be a massive WR 
star rather than a GSPN, given its relatively high apparent brightness, V ~ 11.6 (Massey 1984), 
and the huge interstellar reddening to Ml-67 {Eb-v ~ 0.9-1.5; Solf & Garsenty 1982; Esteban 
et al. 1991; Growther et al. 1995b). This is compatible with My ~ —6.0, which is typical for 
WN8 stars. Photo-ionization modelling of Ml-67 has been made by Growther et al. (1999), who 
were able to trace the central star H-Lyman continuum flux distribution, which is normally not 
directly detectable. The main outcome of this work was to establish the crucial importance of a 
line-blanketed energy distribution for proper modelling of the WN8 atmosphere. 
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As an E-type (ejected material) nebula (see Cliu, Treffers & Kwitter 1983; Smith 1995 and 
references therein), Ml-67 constitutes an object in its earliest phase of wind interaction, when the 
ejecta are least affected by collision and mixing with the interstellar medium, and thus best reflect 
the initial conditions of the ejected winds. This is supported by detailed abundance analysis of the 
nebula: N/0 « 2.95 (Esteban et al. 1991), rather than N/0 « 0.07 for the Galactic ISM (Shaver 
et al. 1983). Note that the E-type Galactic CN RCW58 (central star: WR40; spectral type: 
WN8), is also obviously enriched and must contain stellar ejecta; however this nebula already 
shows a circumstellar ring which suggests a more pronounced interaction with its ambient medium. 
For the same reason, we do not consider the multiple ring (Marston 1995) CN surrounding the 
WN8 star WR16: the two outer rings are likely E-type ring nebulae, but their dynamical ages 
are relatively long (some 10^ yrs and 7 x 10^ yrs). In addition, the inner ring is likely a W-type 
(wind-blown bubble; see Chu et al. 1983 and references therein) CN. 

Ground-based observations {e.g. Chu & Treffers 1981; Sirianni et al. 1998) exhibit a nearly 
spherical, possibly bipolar, patchy structure made up of many bright unresolved knots and 
filaments. The basic properties of Ml-67 are: diameter 110-120" {i.e. 2.4-2.6 pc at a distance of 
4.5 kpc; Grosdidier et al. 1998), mean expansion velocity 40-45 km s^"*^ (Sirianni et al. 1998), 
implying a dynamical age of the order of a few 10^ yrs. 

In the framework of our search for the direct influence of a clumpy, turbulent-like stellar wind on 
the surrounding nebula it ejects and more generally on the interstellar medium (Moffat 1998), 
we have obtained Hubble Space Telescope (HST) Wide Field Planetary Camera Two (WFPC2) 
images of the nebula Ml-67. Figure 1 is an enlargement of the picture reproduced in Grosdidier et 
al. (1998). They reported no overall global shell structure to the nebula (see also Grosdidier et al. 
1999) from the projected radial distribution of the azimuthally-averaged Ha surface brightness. 
Rather, the radial profile centered on WR124 agrees with a decreasing, wind-density power-law 
distribution with increasing distance from the central star. Grosdidier et al. (1998) understood 
the absence of a well defined shell as evidence for an extremely small age of the nebula: the WR 
phase may have just turned on, making the interaction zone {i.e. the nebula) of the present fast 
wind with the previous, slower wind still undetectable (which is in line with the presence of H 
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emission lines in the spectrum of WR124; Crowther et al. 1995a). 

However, these deep HST/WFPC2 images in Ha of the nebula also suggested that Ml-67 may 
likely be the imprint of a previous, luminous blue variable (LBV) slow wind ejected from the WR 
central star WR124, which is supported by other arguments (Sirianni et al. 1998, and references 
therein; Massey & Johnson 1998) as opposed to a classical red supergiant very slow wind (Smith 
1996) blown and compressed by the present WR wind. Note that ubiquitous reversed bow-shock 
like structures are detected throughout the nebula (the more prominent ones being detected at the 
outer edge of Ml-67, PA ~ 10°-20° and 130°-160°; see Figure 2). Generally, the whole periphery 
of the nebula appears made up of numerous small (1-5") reversed bow-shock like structures. 
These features can be interpreted as the result of the thin-shell instability of the highly variable 
previous LBV stellar wind interacting with itself {e.g. Stone et al. 1995). In addition, some dense, 
persisting clumps have possibly been ejected directly from the LBV stellar surface (Grosdidier et 
al. 1998). 

3. Observations and Data Reduction 

3.1. HST Imagery 

Here, we expand on the description given in Grosdidier et al. (1998). Four WFPC2 images 
of Ml-67 with a total combined exposure of 10,000 seconds were taken in March 1997 with the 
HST through the narrowband F656N Hq filter (Biretta et al. 1996). They were combined (with 
the task GCOMBINE in STSDAS/IRAF^) to produce a single deep Ha image cleaned from 
bad/hot pixels and cosmic rays. Before the combination of the four F656N Ha images, we applied 
proper shifts to them in order to achieve alignments to within ±0.01-0.03 pixels. The shifts were 
determined using at least 15 (Planetary Camera; PC) and 25-30 (Wide Field Cameras 1, 2 and 3; 



^IRAF is distributed by the National Optical Astronomy Observatories, operated by the 
Association of Universities for Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 
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WFs) bright stars and calculating the relative displacements between each exposure through the 
estimation of the star centroids. Note that rather than applying one shift to an image to match 
the other, we applied half-displacements to both images. This reduces the differences in artificial 
broadening of the stellar profiles between the two images. Finally, we achieved full widths at half 
maximum spanning 2.4-2.5 pixels (~ 0.11") for the PC, and 1.7-2.0 pixels (0.17-0.2") for the 
WFs, in the Ha band. 

The four narrowband F656N Ha images were straddled in wavelength by broadband V and R 
images: four WFPC2 images of the same field taken through each of the broadband F675W R and 
F555W V filters (Biretta et al. 1996) were separately combined in the same way to produce similar 
high signal-to-noise ratio images of 'continuum' starlight relatively close to the wavelength of Ha. 
The V and R images enabled us to interpolate the stellar light to the wavelength corresponding to 
Ha. After proper flux scaling and spatial shifting (also within ±0.01-0.03 pixels), the interpolated 
broad-band image^'' was carefully subtracted from the Ha image in order to produce a deep 
continuum-subtracted Ha image with the fleld stars removed. Note that we achieved full widths 
at half maximum spanning 2.4-2.5 pixels for the PC, and 1.7-2.0 pixels for the WFs, both in 
the R and V bands. These values are similar to those obtained for the combined Ha image and 
permited an optimum continuum starlight subtraction. 

The scaling in intensity between the interpolated broad-band image and the combined Ha image 
was made by measuring the flux of at least 20 bright stars in each chip, and minimizing the 
residuals in the continuum-subtracted Ha image. This procedure works quite well for Ha line-free 
stars, which are the majority. For WR124 we applied a different flux scaling because this WN8 
star emits a copious amount of Ha+HenA6560 light. It is worthy of note that we have obtained 
all images within a contiguous span of time (6 HST orbits). This reduced problems with potential 
variable stars in the subtraction of starlight from the original Ha image. Note that the images 
were originally pipeline processed before being released but required a full recalibration because of 
severe changes in the dark flies. 
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The V-band component being almost negligible. 
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3.2. Fabry-Perot Interferometry 

Nebular kinematics can be studied very effectively by observations at optical wavelengths 
with the use of a scanning Fabry-Perot interferometer in combination with a two-dimensional 
detector. This technique allows one to probe the velocity field with a wider field of view and 
higher throughput than with a normal long-slit, high-dispersion spectrograph. 
With the servo-stabilized Fabry-Perot (FP) interferometer of the Universite Laval (Joncas & Roy 
1984), we obtained complementary FP Hq data of Ml-67 using MOS/SIS" (Le Fevre et al. 1994) 
at the Canada- France-Hawaii Telescope (CFHT), in 1996 August 30. The instrument was mounted 
at the Cassegrain //8 focus and covered a field of about 200" on a side. An Ha narrow-band filter 
{5\ = 8.8 A) centered at 6575 A was used as a pre-monochromator. The combination of tilting 
and air temperature blue-shifted the central wavelength of the filter to match the red-shifted Ha 
line originating in Ml-67 (mean peculiar velocity: about -1-137 km s~^; Sirianni et al. 1998). The 
FP interferometer has a finesse F = 30 and a free spectral range of 392 km s^^ at Ha, where 
the interference order is p = 765. Since both the interferometer and filter were tilted, we avoided 
ghost images (Georgehn 1970). We used the 2048 x 2048 pixel STIS2 CCD (gain = 4 e" ADU~\ 
readout noise = 8 e~, pixel width = 21 ^) with a spatial sampling of 0.30" per pixel. We employed 
the MOS/SIS tip/tilt image stabilizer in order to improve image quality and achieved effective 
seeing spanning the range FWHM ~ 0.6-0.7". The useful data cube was 1500 x 1500 pixels in 
the spatial domain, and 66 channels (2.2 x F, in order to satisfy the Nyquist criterium) in the 
wavelength domain, giving a velocity sampling of 5.9 km s^^ per channel. The integration time 
was 200 sec channel"^. A neon calibration lamp was used to obtain rest frame interferograms. 
Finally, the data processing was done using the IRAF and ADHOC^^ (see Amram et al. 1995, 
Blais-Ouellette et al. 1999) packages. 



^^http://www. cfht.hawaii.edu/Instruments/Spectroscopy/SIS/ 
^^http://www-obs. cnrs-mrs.fr/ADHOC/adhoc. html 
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4. The Turbulent Status of Ml-67 
4.1. HST Imagery: Wavelets for Structure Function Analysis 

4.1.1. Method 

We favor the use of wavelets to perform structure function calculations because this method 
is much more robust than the classical study of the scaling behavior via the power spectrum 
estimation as a function of the modulus of the wavevector k. Indeed, the latter method is sensitive 
to biases in the data, e.g. smooth polynomial trends, that can easily be shown to drastically affect 
the power spectrum scaling law. In addition, this method is appropriate when the fractal signal 
under consideration is not associated to a unique roughness exponent H (Schmittbuhl et al. 1995). 
Using 2D wavelets, we have analyzed our deep continuum-subtracted Ha image of Ml-67 to isolate 
stochastic structures of different characteristic size. Hence we effectively looked for scaling laws, 
as has been clearly demonstrated in other astrophysical contexts {e.g. clumping in CO: Gill & 
Henriksen 1990; clumping in WR winds: Lepine 1994). In what follows, the deep Ha image will 
be considered as a two-dimensional function noted /. 

Following Muzy et al. (1993), we have explicitly used wavelets to perform a structure function 
analysis. We have considered the field of 'increments' of / over a projected spatial distance r at 
the pixel {x,y), 

A/(r; (x, y)) = f{{x, y) + R/2) - f{{x, y) - R/2), (1) 

by the continuous wavelet transform of / at the scale r, T,j,[/](r; (x,y)), 

A/(r; (x, y)) « n[f]{r; {x, y)) = r'^ x [/ * ^i>r]{r; (x, y)), (2) 

where R is a two-dimensional increment vector (on the plane of the sky) of length r, * denotes 
a convolution and ^r is the analyzing wavelet at the spatial scale r (Muzy et al. 1993). In our 
study, the so-called 'mother-wavelet', ^, is the 2D 'Mexican hat' (formally the second spatial 
derivative of e-(^'+S'')/2: (2 - x^ - y2)e-(^'+j/')/2) ^nd ^^(x, y) = ^(f , f ) (Farge 1992). It is worth 
noting that the Mexican hat wavelet is orthogonal to polynomials of orders and 1. Therefore, 
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the wavelet transform is insensitive to linear trends, hence linear non-stationary components, in 
the signal. 

The statistical moments of |A/(r; {x,y))\ (the two-point correlation known as 'structure function 
of order p') were estimated through the statistical moments of |T>j,[/](r; {x,y))\ by averaging over 
the location: 

{\Afir)n ^ C [ \T^[f]ir;ix,yW dxdy, (3) 



where C is an arbitrary constant. Note that Muzy et al. (1993) applied this formalism to ID 
signals. Here, we generalize their approach to 2D signals. In section 4.1.2 we will perform tests on 
artificial fractal 2D signals, and show that this procedure is correct. Since this approach does not 
depend on the analyzing wavelet (Muzy et al. 1993), one is in a position to test for scaling laws of 
the form: (|A/(r)|P) ~ r'^(p). 

Note that (^(0) = by definition and C,{p) is a smooth, differentiable, monotonically non-decreasing, 
concave function of p (if the signal has absolute bounds), no matter how rough the data are {e.g. 
Marshak et al. 1994). Continuous signals satisfy (|A/(r)|) ~ r, hence C,{p) = p. Processes with 
C{p) Of p are called 'mono-affine' or 'monofractal', whereas processes with variable C{p)/P ^^^ 
called 'multi-affine' or 'multifractal'. 

4.1.2. Tests on Artificial 2D Fractal Signals 

In order to verify the validity of equation (3), it is necessary to test our numerical code on 
artificial 2D signals. Preferably this should be done on 2D fractals in order to check the efficiency 
of our method in finding self-similarity, hence fractality, in the fields of increments. 
The basic method for creating artificial fractals is to consider the spectral density, that is, the 
mean square fiuctuation at any particular frequency and how that varies with frequency. This 
technique, the power spectral method, is based upon the observation that many natural forms 
and signals have a l/u^ frequency power spectrum, i.e. the spectrum falls off as the inverse of 
some power of the frequency i', where the spectral component (3 is related to the fractal dimension 
of the signal: f3 = 2H + 2, with the persistence parameter < H < 1. The latter parameter 
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quantifies the degree of 'roughness' {e.g. H = for a fully rough signal; H = 1 for a fully smooth 
signal) and is related to the fractal dimension, D: H = 3 — D (Moghaddam et al. 1991; Cox & 
Wang 1993). Thus, D will take on values between 2 and 3. 

Figures 3, 4 and 5 show three constructions of artificial fractal signals along with a horizontal cut 
of the images. To build these signals we first generate a random 512 x 512 white-noise signal. We 
then apply a Fourier transform to this image leading to a secondary signal which essentially looks 
like another noise field. We scale the resulting spectra by the desired l/i/^ function and finally 
apply an inverse Fourier transform. Controlling how rough or smooth the surface is has to be 
done by varying the power relationship. From Figure 3 to Figure 5 we have increased the value of 
/3, hence decreasing the roughness. Note that these three images are all based on the same initial 
seed and therefore have the same general shape. 

Figure 6 (left panels) shows our numerical applications of equation (3) for p = \ to the three 
artificial fractal signals. Note that choosing p = 1 in equation (3) and the r~^ scaling in equation 
(2) relates any power-law (|A/(r)|) ~ r^^^> to the persistence parameter H: H = C(l)- Inspection 
of Figure 6 suggests that the conjecture given in equation (3) is very accurate: the derived slopes 
are always very close to the theoretical value of H. 

The right panels of Figure 6 show the hierarchy of exponents C,{p) related to the signals shown in 
Figures 3, 4 and 5. Here the statistical moments of |A/(r; (x,y))| have been estimated for each 
artificial signal, for p = 1-5, in steps of 0.5. The search for scaling laws has been performed as 
before for any individual value of p. Note that our artificial signals are monofractal, therefore 
C,{p) should be strictly proportional to p. This is seen very clearly for H = 0.05. However, for 
H = 0.475 and H = 0.9 we detect a clear departure from this theoretical property for p > 3 and 
p > 1, respectively. This is mainly due to the small amplitude fiuctuations of the last two signals 
(see Figures 4 and 5, lower panels): for large values of p, application of equation (3) may lead to 
very small numbers whenever |Tii,[/](r; {x,y))\ < 1, in combination with dominant extreme events 
which are poorly sampled by definition. 

These numerical experiments demonstrate how our wavelet method is efficient in detecting any 
fractal structure. However, the reliability of the structure function analysis for the highest orders 
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appears very sensitive to the amplitude of the fluctuations. 

4.1.3. Results on the HST/Ha Deep Image of Ml-67 

Figures 7, 8 and 9 show, respectively, the wavelet coefficients of the deep, net Ha image (cf. 
Figure 1) for three particular spatial scales. At the largest scale, 11.3" (Figure 7), the wavelet 
coefficients reveal the large scale structure of Ml-67, where one can easily recognize the bright 
arcs of Figure 1. The number of clumps then increases rapidly towards smaller scales (see Figures 
8 and 9). More strikingly, at the smallest scale close to the resolution limit, 0.3" (Figure 9), the 
density field of Ml-67 appears remarkably structured in apparently chaotically oriented filaments 
everywhere in the nebula. However, many filaments (especially at the smallest scales) seem to be 
slightly radially oriented. Since the field stars have been subtracted off in Figure 1, these filaments 
are strictly related to the fluctuating density field of Ml-67. Note that neither at large scales 
(Figure 7), nor at small scales (Figure 9) is any bipolar outfiow clearly detected. Using the Ha 
HST image. Figure 10 shows the projected angular distribution of the surface brightness averaged 
over 10°-wide sectors centered on WR 124. From the median two possible lobes are detected for 
PA ~ 65°, and PA ~ 185°. Since the lobes are not anti-aligned on the plane of the sky, they are 
unlikely related to a bipolar outflow. Considering the mean, the situation is even worse: many 
subpeaks reveal the highly irregular and patchy structure of Ml-67. Using wider sectors (hence 
losing angular resolution), we arrive at the same conclusions. 

In Figure 11 the quantities (|A/(r)|P) are plotted down to 3 pixels (0.3") of WFPC2, in order 
to secure good sampling of the analyzing wavelet, for p = 1 and 2. (The contribution of the 
instrumental noise has been subtracted off by performing the same analysis on a rectangular 
region outside the nebula; the field stars do not contribute in this calculation since they have been 
removed.) Note the short, marginal inertial range (less than a decade) at the smallest scales from 
about 6.7 X 10^^ to about 1.5 x 10^^ pc which could be attributed to turbulence. More strikingly, 
we detect the presence of a strong scale-break at scales of (L) ~ 5-10" (or 0.11-0.22 pc, adopting 
a distance of 4.5 kpc for Ml-67), corresponding to the characteristic widths of the bright arcs 
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seen in Figure 1. This suggests that the bright arcs may be related to a specific ear her ejection 
event or hydrodynamic instabihty. Looking carefuUy at Figure 11 a second characteristic scale is 
found at (s) ~ 1-2" (or around 3.3 x 10~^ pc), where the second-order structure function is still 
growing but at a significantly lower rate. We interpret these characteristic scales as the presence 
of essentially two types of structures in the signal: bright, small (s)-scale features, and relatively 
less intense, larger (L)-scale features. Indeed, as seen in the second-order structure function which 
emphasizes the most intense events, as the scale r increases (|A/(r)p) first increases sharply 
because the increments from the bottom to the top of the (s)-scale features dominate and there are 
more and more of these. Then (|A/(r)p) becomes saturated with respect to these features: for all 
scales above (s), (|A/(r)p) contributes a relatively fixed number of events to the spatial average. 
Further on we see the larger (L)-scale features at work in a similar way. They allow (|A/(r)p) 
to keep growing but at a slower rate because smaller increments are spawned by these features. 
For scales larger than (L), the second-order structure function again becomes saturated. For 
scales close to the characteristic projected spatial separation between (s)- and (L)-scale features, 
(|A/(r)p) should exhibit a minimum, which is indeed the case for scales around (d) ~ 20", i.e. 
0.4 pc (see Figure 11 and Marshak et al. 1997 for a similar behavior in another physical context). 
Typical voids of width about 20" are compatible with the appearance of Ml-67 in Figure 1. 
Since the structure function of order 1 scales as (|A/(r)|) ~ r^'"^^ {H ~ 0.76) at the smallest scales, 
we infer the fractal dimension of Ml-67 to be about D = 3 — H ^ 2.2-2.3. Thus, Ml-67 appears 
neither as a smooth, nor space-filling structure. For many different nebular objects, we can find 
estimates of their related fractal dimension in the literature. Generally these fractal dimensions 
were derived through box-counting techniques on the perimeter of the nebula, hence leading to 
1 ^ Dp < 2, the majority of the values spanning the range 1.2-1.4 (Bazell & Desert 1988; Blacher 
& Perdang 1990; Falgarone et al. 1991, and references therein). Experiments have shown that 
in the case of isotropic turbulence/fractality, the fractal dimension of the perimeter Dp and the 
fractal dimension D of the two-dimensional projected nebula satisfy the relation: Dp = D — 1 
(Sreenivasan & Meneveau 1986; Beech 1992). Therefore, assuming that the fractal dimension Dp 
of Ml-67 does not vary with the angle of projection (i.e. assuming isotropy of the turbulent status 
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of the nebula), we infer Z)p(Ml-67) ~ 1.2-1.3, which is in good agreement with the measurements 

determined in the above studies and suggests that this parameter takes on universal values. 

In Figure 12 we show the hierarchy of exponents ((p) as a function of p. These exponents have 

been obtained through fitting of the statistical moments over 2.2 decades in scale, including the 

scale break. Note that Cip)/p is no longer constant when p exceeds 2-2.5. We interpret this fact as 

evidence for possible intermittency in the data, i.e. sporadic and intense 'bursts' of high frequency 

activity. 

4.2. CFHT Fabry-Perot Data: Structure Function Analysis 

4-. 2.1. General Results on the Velocity Field 

Figure 13 shows the channels covering the spread in radial velocity for the nebula Ml-67. 
From these data Ml-67 appears as a strongly distorted, thick spherical shell seen almost exactly 
along its direction of rapid spatial motion in the ISM (Moffat et al. 1982; Sirianni et al. 1998). 
The radial velocity of the center of expansion is ~ 137 km s^"*^ (Sirianni et al. 1998). Note that 
towards the star, matter is seen at many different velocities. This is due to the broad Ha emission 
line originating in WR124 (wind terminal velocity: ~ 710 km s~^; Hamann et al. 1993; Crowther 
et al. 1995a). Figure 13 also shows that the velocity dispersion inside Ml-67 is quite large: note 
that a simple thick shell by itself would not explain the multiple radial velocities along the line of 
sight. This velocity dispersion leads one to consider Ml-67 as a thick accelerating shell. On the 
whole, instead of appearing as a nice hollow-type shell projected on the sky, we probably see the 
cap of the bowshock nearly straight on from behind, the far side being greatly intensity-enhanced 
(by a factor of 8 or more) compared to the near side, probably as a result of raming with the ISM. 
This was already claimed by Solf & Carsenty (1982). 

Despite the presence of a 'protective' wind-blown cavity of radius several 10 pc, expected to 
have been produced by the original 0-star progenitor, the current nebula of radius ~ 1 pc will 
nevertheless be able to ram the ISM directly. This is because of the rapid motion of the central 
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star (c. 200 km s^^ ~ 2 x 10^^ pc yr^^) combined with the braking of the cavity expansion in the 
ISM, ahowing the nebula to reach the edge of the cavity and enter in contact with the ISM. 
The processing of ah the interferograms enabled the measurement of the velocity centroids covering 
most of the surface of the nebula. Because of the low signal-to-noise ratio of the approaching, near 
side of Ml-67, this has been done only for the far, 'red' side of the nebula, i.e. for velocities greater 
than 137 km s^^. Figure 14 shows the velocity map of Ml-67 for the red component. From this 
figure, Ml-67 exhibits a distorted velocity field reminiscent of a thick expanding shell as suggested 
by Figure 13. In addition, the velocity field does not show any clear bipolar outflow, contrary to 
the previous claim by Sirianni et al. (1998). We argue that the Sirianni et al. conclusions were 
based on an insufficient amount of data due to poor spatial sampling of the velocity field. Note 
that sudden, huge changes in the transparency conditions of the sky occurred during the scanning 
at the end of the night. This explains the absence of measurements in the extreme south part of 
Ml-67. 

Figure 15 shows the distribution of the 103,287 measured velocity points of the red component. 
The mean LSR velocity is 176.2 km s~^ with a standard deviation of 15.3 km s^^. Given the 
systemic heliocentric radial motion of Ml-67 (~ 137 km s^^) and assuming a symmetric expansion 
of the nebula, with a velocity of expansion of about 46 km s^^ (Sirianni et al. 1998), the maximum 
expected radial velocity should be ~ 183 km s~^. However, we note numerous points above this 
value (up to about 225 km s~^) showing the high velocity dispersion within Ml-67. Note that 
Figure 15 shows velocities well below 137 km s^^, down to about 120 km s^^, i.e. velocities 
apparently related to the blue component. This is mainly due to velocity perturbations of the shell 
near its periphery, and the resulting difficulty in splitting the two components. Moreover, note the 
presence of a small bump near 120-140 km s^^. Since this value is very close to the radial velocity 
of the center of expansion (137 km s^^), here we are mainly concerned with regions of Ml-67 
located near the periphery of the nebula. Such a bump is evidence for velocity perturbations likely 
due to the interaction with the interstellar medium. On the whole, the irregular nature of the 
velocity field is possibly related to either large variations in the density distribution of the ambient 
ISM, or large variations in the central star mass-loss history. 
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4-2.2. Structure Function Analysis of the Velocity Field 

Because any valid statistical analysis of velocity fluctuations has to be done on a (at least 
local) stationary velocity field, one has to subtract global, systematic, nonstationary trends from 
it. Rather than fitting a two-dimensional polynomial to the velocity data in order to remove the 
expansion pattern of the distorted shell (which would generally lead to additional new systematic 
trends), we have convolved the centroid velocity map with Zurflueh filters (Zurflueh 1967; Miesch 
& Bally 1994) of different frequency response widths, and removed the smoothed map from the 
original data. The Zurflueh filter is a symmetric, low-pass filter which nicely approximates a 
step function in the frequency domain and moderates the introduction of artificial high-frequency 
patterns in the convolved data. This systematic trend correction gave us the residual fluctuating 
velocity component on which we have performed the structure function analysis. Figure 16 shows 
the computed second-order discrete structure function of the velocity field of Ml-67: 

where X and r are two-dimensional vectors on the plane of the sky, the summation being 
done on all the pairs separated by r, N(r) (Scalo 1984; Miesch & Bally 1994). In Figure 16, the 
second-order velocity structure function shows no clear, inertial range, as well at small as at larger 
spatial separations. The absence of a well-defined inertial range appears very likely because it is 
independent of the width of the applied Zurflueh filter. However, with a distance of 4.5 kpc to 
Ml-67 and a pixel size of 0.3", we detect a possible, marginally inertial regime occuring in the 
short range 16-50 pixels, i.e. 0.1-0.3 pc. A linear least-squares' fit of the slope in the range where 
there is a possible correlation of the results has been done using a power-law and gives a slope 
of ("(2) ~ 0.3 (i.e. the related distribution of turbulent kinetic energy with spatial scale would 
be described by an energy spectrum, E{k), scaling as k^^'^ with the wave number k). Such a 
value does not meet the Kolmogorov type description for which one would obtain ("(2) = 2/3 (i.e. 
E{k) ~ k^^'"^). On the other hand, the 16-50 pixel possible inertial range is so short that it is 
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unlikely significant or real. 

At smaller scales, down to the achieved spatial resolution, the 2-16 pixel range (i.e. 0.01-0.1 pc) 
seems to be related with an almost stationary, scale-independent field of increments. Indeed, in 
this range S{r) is found to scale trivially: ((2) ~ 0. Note we are not able to confirm the inertial 
range detected from the density field (see section 4.1.3) because this inertial range ends at about 
0.3", which is well below the resolution of the CFHT data. Similar results are obtained when 
estimating the second-order structure function of the centroid velocity field using (linear-trends 
free) wavelet coefficients following the method described in section 4.1.1. 

4-2.3. The Distribution of Velocity Increments in Ml-67 — Structure Function Analysis of the 

Regularized Velocity Map 

Because we depend on the spatial coordinates to perform all averaging operations when 
analysing the velocity data, our first task is to determine which aspects of the signal are more 
likely to be stationary. The velocity field of Ml-67 is non-stationary but with presumably locally 
stationary increments, so that in the previous section we turned to structure function analysis. 
Although possibly related to effects of finite spatial resolution (see Marshak et al. 1994, and 
references therein), the small, almost null (at the smallest scales), value of C,{2) derived from the 
velocity field suggests that the field of increments At'(X; r) = v{X) — v{X + r) in equation (4) is 
apparently essentially scale-independent. 

One-point probability density functions (PDF) make no use of the ordering of the data, leaving 
wide open the possibility of correlations between data-values at different points or structures in the 
data. However, the apparent scale independence of the velocity increments relaxes this drawback 
of one-point PDF statistics. One traditional approach is one-point PDF histograms, typically 
searching for gaussianity. 

Figure 17 shows the empirical PDFs of velocity increments as traced by the velocity field wavelet 
coefficients for four different spatial separations. Each PDF is compared with a Levy stable 
distribution (found using the computer program STABLE; Nolan 1997) and the best-ffi Gaussian 



22 



distribution. Recall that Levy stable distributions arise from the generalization of the central 
limit theorem to a wider class of distributions. In the Gaussian case, the renormalized sum 
of independent and identically distributed random variables with finite variance has the same 
probability distribution as the individual terms in the sum. In 1937 Levy showed that the 
Gaussian distribution is but a special case in a wider family of distributions once the constraint 
of finite variance is relaxed (Feller 1971). The resulting distributions are referred to as Levy 
stable distributions, or simply Levy distributions, and are characterized by slowly decaying tails 
and infinite moments. Levy distributions, Lq,, require four parameters to describe: an index of 
stability a in the interval (0,2], a skewness parameter — 1 < /3 < 1, a scale parameter o" > and 
a location parameter ^. The index a determines the rate at which the tails of the distribution 
decrease {La(x) oc a:~^~"); when a = 2, a Gaussian distribution results. When a < 2 the variance 
is infinite. In general, the p-th moment of a Levy stable random variable is finite if and only if 
p < a. The parameters a and /U determine the width and the shift of the peak of the distribution. 
For a positive (resp. negative) /? the distribution is skewed to the right (resp. the left). For /? = 0, 
the distribution is symmetric about fi. In addition, as a approaches 2, (3 loses its effect and the 
distribution approaches a Gaussian regardless of /3 (Nolan 1997). 

Strictly speaking, velocity increments cannot have a Levy stable distribution because this would 
imply nonzero probability density in regions outside the range for which velocity increments 
are physically meaningful. However, these distributions can yield useful approximations over 
a truncated/finite range as demonstrated in other physical contexts (Painter et al. 1995). In 
addition, when we use an infinite variance distribution in order to fit the distribution of velocity 
increments, all that we are saying is that the observed slowly decaying tails result in behavior 
similar to that of the infinite variance model. On the whole, we stress the fact that variance is but 
one measure of spread for a distribution, and is not appropriate for all problems (Painter et al. 
1995). 

As the scale increases. Figure 17 shows us that the distribution approaches a Gaussian distribution, 
as expected for a truncated Levy distribution (Nakao 2000). At small scales, the distribution of 
the velocity increments is characterized by a non-negligible number of extreme events making the 
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PDF of the velocity increments not narrow enough to enable a simple dimensional argument to 
relate quantitatively all moments: {\Av{r)\'f) no longer approximates {\Av{r)\)'f and precludes 
the determination of any actual (multi-)fractal structure. At large scales the gaussianity of the 
PDFs releases this drawback and non-trivial scaling laws may be detected, although likely biased, 
as seen in Figure 16. Note that the Levy fits in Figure 17 are nearly symmetric {(3 ~ 0) at small 
scales. However, as the scale increases, P increases mainly because the wavelet coefficients take 
into account more and more non-linear trends in the signal. 

Because of the high spatial resolution of our Fabry-Perot observations and the good seeing 
conditions, one expects each pixel to look at a fairly independent piece of the sky. However, very 
low flux regions could see their spectrum polluted by the psf wings of a very bright neighbor. Such 
combination is indeed found in our data and leads to the dominant, larger velocity increments, in 
line with the (s)-scale intense features detected in the density field wavelet analysis (see section 
4.1.3). As a consequence, the velocity centroid measurements are sparsely corrupted resulting 
in fat-tail PDFs of velocity increments. Note that, although there is a regular grid of pixels, 
the emission is distributed irregularly according to the density distribution. Should the velocity 
signal be coming primarily from the brightest regions then the scaling signal could be additionally 
confused by the irregular grid signal (see Rauzy et al. 1995). However, our FP, high throughput 
velocity map is not found to suffer from a significant interference signal from a non-uniform grid. 
We therefore applied a 1-pixel radius ring median filter to the velocity map. The ring median 
filter is defined as a median filter which assigns weights only to selected pixels in an annulus. Its 
advantage is that it has a sharply-defined scale length; that is, all objects with a scale-size less 
than the radius of the ring are filtered and replaced by the local background level. It provides a 
simple method to remove noise-like deviations (independent of morphology) from a digital image, 
leaving behind the large-scale structures and overall gradients (Seeker 1995). We expect each 
velocity point of the regularized velocity map to be an accurate measure of the underlying gas 
kinematics, averaged on its pixel seeing disk. 

Figure 18 shows the computed structure functions of the velocity field of Ml-67 after small-scale 
regularization using a 1-pixel radius ring median filter. The scaling range [~ 1.4", 10"] where the 
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exponents are defined is clearly marked. The break which defines the upper limit of ~ 10" is not 
robust with respect to adding to the average more data points, and simply reflects a violation of 
the ergodicity assumption. Below the 1.4" scale, the slope is slightly steeper, as a consequence of 
the ring median filtering which tends to smooth the signal at the smallest scales. Indeed, using 
a 2-pixel ring annulus simply results in moving the lower limit for scaling towards larger scales. 
Therefore, the cleaned velocity map is found to scale over at least one order of magnitude in scale, 
with C(l) ~ 0.48-0.49 and C(2) ^ 0.90-0.91. The latter value is close to 1.0 {E{k) ~ k''^-^), 
suggesting the dynamics is essentially shock dominated (for which E{k) ~ fc~^). Figure 19 shows 
the function (^ as a function of p, up to 8.0, in steps of 0.1. Clearly C{p)/p is non-linear and 
suggests the multi-afhnity of the velocity field. 

Multi-affine modelling could be done in the context of 'Universal Multifractals' (Schertzer & 
Lovejoy 1987) where multifractals are the generic result of multiplicative cascades. A continuous- 
scale limit of such processes leads to the family of log-infinitely divisible distributions which have 
a Gaussian or Levy stable generator (Schertzer & Lovejoy 1987, 1997; Lovejoy &: Schertzer 1990). 
For universal multifractals we have: 

Cip)=pC{l)--^{p--p) (a /I), (5) 

a — I 

dp) = pCil) - CipHp) (a = l), (6) 

where Ci < 2 is an intermittency parameter and < a < 2 is the Levy tail index. After the 
existence of a scaling regime is established, one can test for universal multifractal behavior. An 
efficient technique for that purpose is the 'Double Trace Moment' (DTM) method (Lavallee 1991). 
A new function K{q,ri) is first defined as K{qr]) — qK{rj), with K{q) = gC(l) ~ CiQ)- ^^ the case of 
universal multifractals, the ry-dependence in K{q,ri) factorizes as: K{q,r]) = rj'^K^q). The DTM 
method allows the computation of K(q,r]), and assuming universality, the Levy index a can be 
estimated by fixing q and varying rj. With a good estimate of a in hand, an ordinary least-squares 
estimate of Ci is easy to obtain. In our case, a standard two-parameter nonlinear regression does 
very fine over almost two orders of magnitude in rj, with a ~ 1.91 and Ci ~ 0.04 it 0.01 (see 
Figure 20). The related C{p) universal multifractal fit is shown in Figure 19 (solid curve). The 
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agreement between the empirical ({p) and the multifractal fit is very good for p below ~ 7. Note 
that the finite size of the sample implies that sufficiently high order moments are dominated by 
the largest values in the field and therefore underestimate the true ensemble moments. Beyond 
uiayi{qrj,r]) = {2/CiY''^ ^ "^-"^ equation (5) is no longer true, and C,{p) becomes linear (see Figure 
19). This is the so-called first order multifractal transition (Schertzer & Lovejoy 1992) which is 
also seen in Figure 20 where the empirical K{q^ rjYs are no longer fitted by the linear regressions. 

5. Discussion and Conclusion 

Using the unique, high-resolution, wide-field imaging capabilities of HST, we have tested 
quantitatively, via scaling laws, for compressible turbulence in the distribution of density clumps in 
the ejected nebula Ml-67. These data have been straddled by complementary CFHT Fabry-Perot 
Ha data in order to perform the same test on the velocity field of Ml-67. What is generally 
understood as turbulence is the difference resulting from the subtraction of the instrumental 
and thermal broadening from the observed emission line widths. In other words, turbulence is a 
quantity related to macroscopic, chaotic motions of the gas, and their related density fluctuations. 
At small scales, the density field of Ml-67 appears remarkably structured in essentially chaotically 
oriented filaments everywhere in the nebula. However, many filaments (especially at the smallest 
scales) seem to be slightly radially oriented. This is not surprising, since the most likely force 
which drives the nebula and the ensuing turbulence is directed outwards in all directions from the 
central star. The presence of filaments agrees with earlier numerical simulations of turbulent flows 
{e.g. Sreenivasan & Antonia 1997; Vazquez-Semadeni 1999, and references therein) or observations 
{e.g. Miville-Deschenes et al. 1999, and references therein) which have revealed how turbulence 
is characterized by fllamentary vortical structures and surrounding dissipative sheets. It is worth 
noting that the empirical PDFs in Figure 17 are statistically good estimates of the true PDFs in 
the central parts of the distributions. Clearly, these PDFs are not gaussian at the smallest scales; 
rather, they are best described by Levy stable distributions which are characterized by extended 
tails. In addition, as the spatial scale increases, the Levy index a is found to increase, the gaussian 
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case (q = 2) being attained for scales above about 4". This behavior also suggests the presence of 
vorticity in the flow and is likely related to the phenomenon of intermittency (Lis et al. 1998, and 
references therein). 

The velocity field of Ml-67 (after correction for bad velocity points related to the largest velocity 
increments) shows a clear inertial range reminiscent of turbulence. Because the ISM and the 
ejected nebula are compressible and inhomogeneous, the gas flows being in addition highly 
supersonic, we expect shock wave propagations within Ml-67. External forces due to stellar winds 
(or gravitation in another astrophysical context like Giant Molecular Clouds: Gill & Henriksen 
1990) play an important role in feeding the ISM with highly pressurized motions, moving the gas 
at velocities larger than the sound speed, and ultimately leading to supersonic turbulence and its 
related small scale compressions. Interaction between shock waves and turbulent fluctuations are 
expected to distribute the energy dissipation on a cascade characterized by a steeper law ^(1) = 0.5 
(Fleck 1996) than predicted by Kolmogorov's law (C(l) = 1/3). Our value of C(l) = 0.48-0.49 
strongly suggests that the turbulence traced by the velocity field is compressible. In addition, since 
our observations come from a projection on the plane of the sky, we therefore expect a smearing 
of the information along each line of sight. This may lead to additional perturbations in the flows, 
with apparently fluctuating power-laws as a function of projected spatial separation (O'Dell & 
Castaheda 1987, and references therein) precluding any inertial regime detection. Figure 18 tells 
us that this effect is practically negligible, and additionally suggests that the Ha layer is thin 
compared to the projected spatial separations. On the other hand, the value C,{2) = 0.90 is in 
very good agreement with the empirical Larson-type law which predicts Q{2) = 0.8 it 0.1 (Larson 
1981; Miesch &: Bally 1994) and is associated to fields which are essentially shock-dominated. 
Note that C,{2) < 2 x ("(1) results from the slight intermittency of the velocity field, which 
makes it multifractal rather than monofractal. With Ci close to 0, the process is found to be 
nearly homogeneous (the fractal dimension of the set contributing to the mean velocity field is 
2 — Ci ~ 1.96). On the other hand, the degree of multifractality a ~ 1.90-1.92 is close to 2. The 
latter result suggests that we are far from the so-called monofractal /3-model (a = 0); rather, the 
turbulence is best described by the log-normal model {a = 2) for which high values of the field 
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dominate more than for smaller values of a. The estimation of a and Ci for other HII regions 
appears now essential in order to test whether the multifractal parameters found in the case of 
Ml-67 take on universal values. For comparison with other turbulence studies, one may derive 
from Ci the standard intermittency parameter ^ which is the autocorrelation exponent of the 
dissipation field: ^ = K{2, 1). For the /3-model (q = 0), n = Ci, whereas for the lognormal model 
{a = 2), n = 2 X Ci. In our case, with a « 1.91, we get ^ ^ 1.90 x Ci ^ 0.076. 
Recall that turbulent velocity measurements in wind tunnel experiments lead to a Kolmogorov-type 
(i.e. incompressible; C(2) = 2/3) turbulence with a ^ 1.3±0.1, Ci « 0.25±0.05 and fi » 0.35±0.1 
(Schmitt et al. 1992). The higher degree of multifractality and the smaller intermittency found 
for the velocity field of Ml-67 is likely a consequence of the compressible nature of the tracing gas. 
This will be discussed in a forthcoming paper. 

It is worthy of note that large variations may exist in the history of WR124's mass-loss and 
the distribution of the ambient ISM. This precludes 1) any unique interpretation of the absence 
of a well-defined inertial range in the density field; and 2) any unambiguous derivation of the 
stellar mass-loss history because density fiuctuations in the absence of self-gravitational effects 
are believed to be transient only in the context of supersonic turbulence (Vazquez-Semadeni 
1994). However, note that the marginally inertial range found from the density field is likely a 
consequence of the rapid regeneration of density structures triggered by the turbulent nature of 
the velocity field. Finally, given the extreme perturbations of the velocity and density fields in 
Ml-67, it is virtually impossible to measure any systematic impact of the WR (or LBV) wind on 
the nebular structure. For whatever reason, any systematic effects seem to have been randomized 
such that they have become fully masked. On the other hand, the origin of the (s)-scale features 
found from the wavelet analysis of the HST Ha image still has to be investigated in detail. These 
features could be related to either clumps of stellar origin, or smoothed out versions of the bright 
'bullets' reported in Grosdidier et al. (1998). On the other hand, the irregular nature of the 
velocity field is likely due to either large variations in the density distribution of the ambient ISM, 
or large variations in the central star mass-loss history. 
On the whole, Ml-67 exhibits the signatures of compressible, intermittent turbulence, the driver 
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being ultimately the stellar wind originating in WR124, as opposed to gravity/magnetic fields in 
the ISM. In addition, from the structure function of order 1 of the density field (which scales as 
(|A/(r)|) ~ 7-0-76 ^^ ^]^g smallest scales), we infer the fractal dimension of Ml-67 to be D ~ 2.2-2.3. 
This result is very close to the values often quoted for molecular clouds and indeed is essentially 
the result found in the wavelet analysis of Gill & Henriksen (1990). Therefore, although different 
drivers may be related to the generation of supersonic, compressible turbulence, it appears that 
the resulting phenomenology essentially remains the same. 
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Fig. 1. — Deep WFPC2/Ha log-scale image of Ml-67 from 4 combined exposures totalling nearly 
3 hours. The field stars have been subtracted out (enlargement of Figure 1 from Grosdidier et al. 
1998). 
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Fig. 2. — Enlargement of the faint outer north-eastern region {PA « 10°-30°) of the nebula (upper 
panel) where at least 9 'reversed' bow-shocks are visible (see lower panel). Such features are seen 
throughout the nebula although better detected at the periphery. A few residuals of subtracted 
field stars also appear. The chosen contrast emphasizes the dimmer features. 
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Fig. 3. — Artificial fractal signal with H = 0.05 (upper panel), along with a horizontal cut in the 
middle of the image (lower panel) . 
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Fig. 4. — Artificial fractal signal with H = 0.475 (upper panel), along with a horizontal cut in the 
middle of the image (lower panel) . 
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Fig. 5. — Artificial fractal signal with H = 0.90 (upper panel), along with a horizontal cut in the 
middle of the image (lower panel) . 
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Fig. 6. — Left panels: Persistence parameter retrieval for H = 0.05, H = 0.475 and H = 0.9 (cf. 
Figures 3, 4 and 5). Right panels: Structure function analysis of the signals shown in Figures 3, 4 
and 5. 
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Fig. 7. — Wavelet coefficients of Ml-67 at tlie scale 11.3". The orientation of the image is the same 
as of Figure 1. 
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Fig. 8. — Wavelet coefficients of Ml-67 at the scale 1.8". The orientation of the image is the same 
as of Figure 1. 
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Fig. 9. — Wavelet coefficients of Ml-67 at the scale 0.3". The orientation of the image is the same 
as of Figure 1. 
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Fig. 10. — The projected angular distribution of the surface brightness averaged over 10°-wide 
sectors centered on WR 124. 
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Fig. 11. — Structure function analysis of Ml-67. The quantities (|A/(r)|^) are plotted down to 3 
pixels of WFPC2 (0.3"), for p = 1 and 2. Three characteristic scales are also indicated (see text). 
One arcsecond corresponds to about 2.2 x 10~^ pc assuming a distance of 4.5 kpc for Ml-67. 
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Fig. 12. — Structure function analysis of Ml-67. The corresponding ("(p) function demonstrates 
the multi-affinity of Ml-67. 
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Fig. 13. — Maps of the CFHT/SIS-Ha intensity in Ml-67 for the heUocentric radial velocities (in 
km s^^) indicated in the upper right-hand corners. Everywhere in the nebula except near the edges, 
a splitting of the Ha line into at least two components is detected; the high-velocity component is 
generally brighter (by a factor 8 or more) than the low-velocity component. The central star shows 
up clearly at the same position in each box. 
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Fig. 14. — Smoothed velocity map of Ml-67 for the hehocentric radial velocities of the 'red' 
component (in km s~^). North is up, East to the left. Note the missing chunk of emission in 
the bottom part due to sudden arrival of clouds near the end of the FP scan. 
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Fig. 15. — The distribution of the 103,287 measured velocity centroids of the red component. 
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Fig. 16. — The second-order velocity structure function of Ml-67 where LAG is the separation 
between velocity pairs. One pixel corresponds to ~ 6.5 x 10~^ pc assuming a distance of 4.5 kpc 
for Ml-67. 
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Fig. 17. — Empirical PDFs of velocity increments as traced by the velocity field wavelet coefficients 
for four different spatial scale separations (squares). Each PDF is compared with a Levy stable 
distribution (solid curve) and the best-fit Gaussian distribution (dotted curve). The histograms 
were constructed from about 27300 individual measurements at each scale. 
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Fig. 18. — Velocity structure function analysis of Ml-67 after small-scale regularization using a 1- 
pixel radius ring median filter. The quantities {\Av{r)\P) are plotted down to 3 pixels of CFHT/SIS 
(0.9"), for p = 0.5-2.5, in steps of 0.5. One arcsecond corresponds to about 2.2 x 10^^ pc assuming 
a distance of 4.5 kpc for Ml-67. 
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Fig. 19. — Velocity structure function analysis of Ml-67 after small-scale regularization using a 
1-pixel radius ring median filter. The corresponding ("(p) function (points) demonstrates the multi- 
affinity of the velocity field of Ml-67. A universal multifractal fit is also shown as a solid curve (see 
text). 
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Fig. 20.— Plots of logio \K{q,r])\ (using q = 0.5, 1.5, 2, 2.5, 3 and 5.5) obtained with a DTM 
analysis of the velocity field of Ml-67 after small-scale regularization. From the slopes and the 
logio ^ = intercepts of the linear region we obtain a = 1.90-1.92, and Ci = 0.04 it 0.01. There 
is a first order multifractal phase transition near uiax(qr], rj) = 7, which is a consequence of the 
finiteness of the analyzed sample. 



